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We perform numerical simulation based on the time-dependent mean-field Gross-Pitaevskii equa- 
tion to understand some aspects of a recent experiment by Donley et al. on the dynamics of 
collapsing and exploding Bose-Einstein condensates of ®^Rb atoms. They manipulated the atomic 
interaction by an external magnetic field via a Feshbach resonance, thus changing the repulsive 
condensate into an attractive one and vice versa. In the actual experiment they changed suddenly 
the scattering length of atomic interaction from positive to a large negative value on a pre-formed 
condensate in an axially symmetric trap. Consequently, the condensate collapses and ejects atoms 
via explosion. We find that the present mean-field analysis can explain some aspects of the dynamics 
of the collapsing and exploding Bose-Einstein condensates. 
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I. INTRODUCTION 

Recent successful detection || of Bose-Einstein 

condensates (BEC) in dilute bosonic atoms employing 
magnetic trap at ultra-low temperature has intensified 
experimental activities on various aspects of the con- 
densate. On the theoretical front, numerical simulation 
based on the time-dependent nonlinear mean-field Gross- 
Pitaevskii (GP) equation Q has provided a satisfactory 
account of some of these experiments ^ 0, § ■ Since 
the detection of BEC for ^Li atoms with attractive inter- 
action, one problem of extreme interest is the dynamical 
study of the formation and decay of BEC for attractive 
atomic interaction 

For attractive interaction the condensate is stable for 
a maximum critical number A^cr of atoms [||. Measure- 
ments for A^cr 1,0 are in reasonable agreement with the 
mean-field analyses for BEC of ^Li in a spherically sym- 
metric trap P, although there is some discrepancy 
for BEC of ^°Rb in an axially symmetric trap ^ . If 
the number of atoms can somehow be increased beyond 
this critical number, due to interatomic attraction the 
condensate collapses emitting atoms until the number of 
atoms is reduced below A^cr and a stable configuration is 
reached. With a supply of atoms from an external source 
the condensate can grow again and thus a series of col- 
lapses can take place, that was observed experimentally 
in the BEC of ^Li with attractive interaction Theo- 
retical mean-field analysis has been able to explain this 
dynamics [||, |, 0, |, §• 

Recently a more challenging experiment has been per- 
formed by Donley et al. at JILA on an attractive 
condensate of ^^Rb atoms [pi)[ in an axially symmetric 
trap, where they manipulated the interatomic interac- 
tion by changing the external magnetic field exploiting a 
nearby Feshbach resonance . In the vicinity of a Fes- 
hbach resonance the atomic scattering length a can be 
varied over a huge range by adjusting the external mag- 
netic field. Consequently, they were able to change sud- 



denly the atomic scattering length by a large amount for 
a BEC of *^Rb atoms They even changed the sign 
of the scattering length, thus transforming a repulsive 
condensate into an attractive one. The original experi- 
ment on attractive ^Li atoms did not use a Feshbach 
resonance, hence the atomic interaction was fixed. This 
restricts the number of atoms in the ^Li BEC to a num- 
ber close to A'cr and the collapse is driven by a stochastic 
process jH p^ . 

In the experiment conducted at JILA, Donley et al. 
changed a stable pre-formed repulsive condensate of *^Rb 
atoms into a highly explosive and collapsing attrac- 
tive condensate and studied the dynamics of collapsing 
and exploding condensates The natural scattering 

length of *^Rb atoms is negative (attractive). By exploit- 
ing the Feshbach resonance they made it positive (repul- 
sive) in the initial state, where the number of atoms, un- 
like in the experiment with ^Li [^, could be arbitrarily 
large. So immediately after the jump in the scattering 
length to a large negative value, one has a highly un- 
stable BEC, where the number of atoms could be much 
larger than A^cr- Donley et al. have provided a quantita- 
tive estimate of the explosion of this BEC by measuring 
the number of atoms remaining in the condensate as a 
function of time until an equilibrium is reached. They 
claim that their experiment reveal many interesting phe- 
nomena that challenge theoretical models. Because the 
phenomenon looks very much like a tiny supernova, or 
exploding star, the researchers dubbed it a "Bosenova" . 
The fundamental physical process underlying the explo- 
sion remains a mystery. 

In this paper we perform a mean-field analysis based 
on the time-dependent GP equation to understand some 
aspects of the above collapse and explosion of the attrac- 
tive condensate of *^Rb atoms in an axially symmetric 
trap. To account for the loss of atoms from the strongly 
attractive condensate we include an absorptive nonlin- 
ear three-body recombination term in the GP equation. 
Three-body recombination leads to the formation of di- 
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atomic molecules with liberation of energy responsible for 
energetic explosion with ejection of matter from the BEC. 
This process could be termed "atomic fusion" in contrast 
to nuclear fusion in stars. The three-body recombination 
rate wc use in numerical simulation is in agreement with 
previous experimental measurement and theoretical 
calculation jl^. The numerical method, we use, for the 
solution of the time-dependent GP equation with an ax- 
ially symmetric trap has appeared elsewhere ^ . 

We find that the present mean-field numerical simula- 
tion provides a fair description of some features of the 
experiment at JILA [l^ ]. 

There have been other theoretical studies based on the 
mean-field GP equation H, 0, |§1 to deal with dynam- 
ical collapse including an absorptive term to account for 
the loss of particles. In most cases the loss mechanism is 
three-body recombination as in the present study. How- 
ever, Duine and Stoof [Q propose that the loss arises 
due to a new elastic process. Instead of attempting a full 
numerical solution of the GP equation with axial symme- 
try, these investigations used various approximations to 
study the time evolution of the condensate or ernployed 
a spherically symmetric trap. Duine and Stoof con- 
sider the full anisotropic dynamics, but use a Gaussian 
approximation for the wave function rather than an exact 
numerical solution. Most of the other studies employed 
a spherically symmetric trap ^. However, the inves- 
tigation of Ref. ^ employed an axially symmetric trap 
to describe some aspects of the experiment at JILA and 
we comment on this work in Sec. IV. In the present in- 
vestigation we consider the complete numerical solution 
of the mean-field GP equation for an axially symmetric 
trap as in the experiment at JILA. It is realized that an 
approximate solution as in the previous studies can not 
explain the dynamics of this experiment . 

In Sec. II we present the theoretical model and the 
numerical method for its solution. In Sec, HI we present 
our results that we compare with the experiment at JILA. 
Finally, in Sees. IV and V we present a brief discussion 
and concluding remarks. 



II. NONLINEAR GROSS-PITAEVSKII 
EQUATION 

A. Theoretical Model Equations 

The time-dependent Bose-Einstein condensate wave 
function \l/(r;T) at position r and time r allowing for 
atomic loss may be described by the following mean-field 
nonlinear GP equation 0, H] 



Here m is the mass and A'' the number of atoms in the 
condensate, g = Ai:h?a/m the strength of interatomic in- 
teraction, with a the atomic scattering length. A positive 
a corresponds to a repulsive interaction and a negative 
a to an attractive interaction. The terms K2 and 
denote two-body dipolar and three-body recombination 
loss-rate coefficients, respectively. There are many ways 
to account for the loss mechanism 0. It is quite im- 
possible to include them all in a self consistent fashion. 
Here we simulate the atom loss via the most important 
quintic three-body term H Jll' contribution 

of the cubic two-body loss term is expected to be 
neghgible ||] compared to the three-body term in the 
present problem of the collapsed condensate with large 
density and will not be considered here. 

The trap potential with cylindrical symmetry may be 
written as V{y) = ^'muj^{r'^ + A^z^) where u is the an- 
gular frequency in the radial direction r and Aw that in 
the axial direction z. We are using the cylindrical co- 
ordinate system r = (r, 0, z) with 9 the azimuthal an- 
gle. The normalization condition of the wave function is 
/dr|*(r;r)|2 = 1. 

In the absence of angular momentum the wave function 
has the form \E'(r; r) = ^/;(r, z; r). Now transforming to di- 
mensionless variables defined by a; = y/2r/l, y = ^/2z/l, 
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where n = Na/l and ^ = AK^/ [a^l^oj). This scaled mean- 
field equation has the correct n dependence of the three- 
body term so that the same equation can be used to 
study the decay rate of different initial and final scatter- 
ing lengths ajni^jal ^^nd ^collapse' respectively, and initial 
number of atoms Nq. In this study the term ATa will be 
used for a description of atom loss in the case of attrac- 
tive interaction, where the scattering length a is negative. 
From theoretical |^ and experimental studies it has 
been found that for negative a, Kj, increases rapidly as 
|a|", where the theoretical study favors n = 2 for 
smaller values of \a\. For larger la |, a much larger rate 
of increase may take place There are theoret- 

ical y, M ml and experimental estimates of 
for ^'Rb,™Na, and ^Li away from Feshbach resonance. 
However, no thorough and systematic study of the varia- 
tion of ATs near a Feshbach resonance has been performed 
p3[ . An accurate representation of the variation of K^^ 
of^^Rb near the Feshbach resonance is beyond the scope 
of this study and here we represent this variation via a 
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quadratic dependence: K3 ~ a^. This makes the param- 
eter ^ above a constant for an experimental set up with 
fixed I and lu and in the present study we use a constant 
t 

The normahzation condition of the wave function be- 
comes 

AAnorm = 2^/ dx dy\^{x,y;t)\^x-^ = 1. (2.4) 

For K3 ~ 0, A/norm — 1, however, in the presence of loss 
K3 > 0, A/norm < 1- The number of remaining atoms N 
in the condensate is given by iV = A^oA/norm, where A^o 
is the initial number. 

The root mean square (rms) sizes Xrms and yrms are 
defined by 

poo poo 

a^rms = -^110^27^ / dx dy\ip{x,y]t)\'^x,{2.5) 

Jo J ~oo 

POO POC 

J/ms = •^noVm27r / dx dy\ip{x,y;t)\'^y'^x~^ (2.6) 

Jq J -oc 

B. Numerical Detail 



using a 



We solve the GP equation (2.3) numerically 
time- iteration method elaborated in Refs. |l8 , 19 , 24| . 
The full GP Hamiltonian is conveniently broken into 
three parts — H^, Hy, and _ff„ — the first contain- 
ing the x-dependent linear terms, the second contain- 
ing the y-dependent linear terms and the third contain- 
ing the nonlinear terms. The GP equations for the first 
two parts are defined on a two-dimensional set of grid 
points Nx X Ny using the Crank-Nicholson discretization 
method. The resultant tridiagonal equations along x and 
y directions are solved alternately by the Gaussian elim- 
ination method along the x and y directions [^. The 
GP equation for the third part do not contain any space 
derivative and is solved essentially exactly. Effectively, 
each time iteration of the GP equation is broken up into 
three parts — using H^, Hy and i?„. For a small time step 
A the error involved in this break-up procedure along 
X and y directions is quadratic in A and hence can be 
neglected. For numerical purpose we discretize the GP 
equation using time step A = 0.001 and space step 0.1 
for both X and y spanning x from to 15 and y from —30 
to 30. This domain of space was sufficient to encompass 
the whole condensate wave function even during and af- 
ter collapse and explosion. The preparation of the initial 
repulsive wave function is now a routine job and was 
done by increasing the nonlinearity n of the GP equation 
(2.3) by 0.0001 in each time step A during time iteration 
start ing with the known harmonic oscillator solution of 
Eq. (IJ) for n = ^ = @. 

It is now appropriate to calculate the parameters of 
the present dimensionless GP equation ( p. 3D correspond- 
ing to the experiment at JILA. We follow the notation 
and nomenclature of Ref. [ p^ . Their radial and axial 



trap frequencies are i^j-adial ~ ^^-^ ^^'^ '^axial ~ ^-^ 
Hz, respectively, leading to A = 0.389. The harmonic 
oscillator length / of ^^Kh ato ms for u; — 2tt x 17.5 Hz 
and m « 79176 MeV is I = y/h/{muj) = 26070 A. One 
unit of time t of Eq. (U) is 1/w or 0.009095 s. They 
prepared a stable ^^Rb condensate of A^o = 16000 atoms 
with scattering length ojnitial ~ ''''^0, Qo = 0.5292 A, 
such that the initial n = 2.274. Then during an inter- 
val of time 0.1 ms the scattering length was ramped to 

~ '^collapse ~ — 30ao such that final n = —9.744. The 
final condensate is strongly attractive and unstable and 
undergoes a sequence of collapse and explosion. 

The initial value of n(= 2.274) was attained after 22740 
time steps. The nonlinearity n is then ramped from 2.274 
to —9.744 in 0.1 ms. As one unit of dimensionless time 
t is 0.009095 s, 0.1 ms corresponds to 11 steps of time 
A. In the present simulation, n is ramped from 2.274 
to —9.744 in the GP equation by equal amount in 11 
steps. The absorptive term ^ was set equal to zero during 
above time iteration. Now the system is prepared for the 
simulation of the collapse and explosion. 

For the simulation of the collapse and explosion the cu- 
bic nonlinear term is maintained constant and a nonzero 
value of ^ is chosen. The time-evolution of the GP equa- 
tion is continued as a function of time t = Tg^Qj^g start- 
ing at 0. The time-evolution is continued using time step 
A = 0.001. After a small experimentation it is found 
that ^ = 2 fits the experiment at JILA satisfactorily. 
Unless otherwise specified, this value of ^ was used in all 
simulations reported in this paper for different ^initial' 
^coUapse' ^nd A^o- 

It is useful to compare this value of £,{= 2) with the 
experimental ||l^ and theoretical |Q estimates of three- 
body loss rate of ^^Rb. For this we recall that K3 = 
S^o?l^u}/A. Under experimental condition of an external 
magnetic field of 250 gauss on ®^Rb the scattering 
length was a ^ — 370ao. Consequently, the present value 
of 2) corresponds to Kj, ~ 9 x 10~^^ cm^/s for a ~ 
— 370ao, which is about two times the experimental rate 
A'3 = (4.24l°;™ ± 0.85) x lO^^s gj^6/g g ^j^^^^^ 

1.3 times the theoretical rate Kj, = 6.7 x 10^^'°^ cm^/s at 
a~ -370ao 0. 



III. NUMERICAL RESULT 



The numerical simulation using Eq. (2^) with a 
nonzero ^ immediately yields the remaining number of 
atoms in the condensate after the jump in scattering 
length. The remaining number of atoms vs. time is plot- 
ted in Fig. 1 for ajnitial = 7ao, flcollapse = -30ao, 
^ = 2, and A^o = 16000 and compared with the exper- 
imental data. In this figure we also plot the result in 
this case for ^ = 3, which leads to a better agreement 
with experiment for this specific case. However, the use 
of ^ = 2 leads to a more satisfactory overall agreement 
with experiment. Except this single curve in Fig. 1 and 
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the plot in Fig. 4 (a) below, which are calculated with 
^ = 3, all results reported in this paper are calculated 
with C = 2. 
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FIG. 1: Number of remaining atoms in the condensate of 
16000 *^Rb atoms after ramping the scattering length from 
"initial ~ "collapse ~ — 6.7ao, and — 30ao, — 263ao in 

0.1 ms as a function of evolution time Taifolve ^^lid 
circle: experiment for ctgoUapse ~ — 30ao [ [l4| ; full line: theory 
= 2); dash-dot line: theory = 3, dcoUapse ~ — 30ao); 
dash ed l ine: average over preliminary, unanalized data using 
Eq. Hi. 

In the experiment at JILA it was observed that the 
strongly attractive condensate after preparation remains 
stable with a constant number of atoms for an interval 
of time ^collapse' *^^ll*5d collapse time. This behavior 
is physically expected. Immediately after the jump in 
scattering length from 7ao to — 30ao, the attractive con- 
densate shrinks in size during ^collapse' central 
density increases to a maximum. Then the absorptive 
three-body term takes full control to initiate the explo- 
sion. Consequently, the number of atoms remains con- 
stant for T^^QYye < '•collapse- present result (full 
line) also shows a similar behavior. However, in this sim- 
ulation the absorptive term is operative from Tgyolve ~ ^ 
and the atom number decreases right from beginning, al- 
beit at a much smaller rate for Tgypj^g < 'collapse' 

Donley et al. repeated their experiment with different 
values of a-^-^-^^, Ocollapse' No For a-^-^-^^ - 
Tag we repeated our calculation with the following val- 
ues of final scattering length: Ocollapse ^ — 263ao and 
— 6.7ao. These results are also plotted in Fig. 1 and 
agree with the unpublished, preliminary unanalyzed data 
p5| . The initial delay ^collapse ™ starting the explosion 



is large for small | ^(jQj^g^pgg | as we see in Fig. 1. Sim- 
ilar effect was observed in the experiment for an initial 
condensate of 6000 atoms as shown in their Fig. 2 p[ . 
After a sequence of collapse and explosion, Donley et al. 
observed a "remnant" condensate of A^remnant atoms at 
large times containing a certain constant fraction of the 
initial Nq atoms. Figure 1 shows such a behavior. 
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FIG. 2: The central part of the dimensionless wave function 
\(f>{x,y)\ = \ip{x,y)/x\ of the condensate on 0.1 x 0.1 grid for 
^ = 2 after the jump in the scattering length of a BEC of 
16000 *^Rb atoms from ajj^j^jai ~ 7ao to ffl^ollapse ~ — 30ao 
at times TgyQ^^g = 0, 3.6 ms, 3.8 ms, and 8 ms. The quantities 
X and y are expressed in units of //\/2, where I = 26070 A. 

The above evolution of the condensate after the jump 
in scattering length to — 30ao from 7ao for Nq = 16000 
can be understood from a study of the wave function 
and we display the central part of the wave function in 
Fig. 2 for t^yqIyq = 0, 3.6, 3.8, and 8 ms. The wave 
function immediately after jump at time ''"gyolve = is 
essentially the same as that before the jump at —0.1 ms. 
There is not enough time for the wave function to get 
modified at TgyQjyg = 0. From Fig. 2 we find that at 
3.6 ms the wave function is only slightly narrower than 
at ms but still smooth and has not yet collapsed suf- 
ficiently. As TqyqIyq increases, the wave function con- 
tracts further and the explosion starts. At 3.8 ms some 
spikes (irregularities) have appeared in the wave function 
showing the beginning of the explosion and loss. From 
the study of the wave functions we find that the explo- 
sion start at Tgypj^g = 'collapse — '^ '^ agreement 
with the experiment at JILA. We also find that at 3.7 
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ms before the loss began the bulk BEC did not con- 
tract dramatically as also observed in the experiment. 
In the numerical simulation for this case we find that at 
"^evolve ~ Oj2;rms = 2.98^m and j/rms — 4.21^m and at 
■^evolve = ^rms = 2.53^m and j/rms = 4.10^m. 

From Fig. 2 we see that at 8 ms the wave function is very 
spiky corresponding to the violent ongoing explosion. 
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Donley et al. fitted the decay m the number of atoms 
during particle loss to a decay constant T^^^f^^y via the 
formula 
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FIG. 3: More decay curves with ^ = 2 for (a) ajnitial ~ 89ao, 
"collapse ^ -15ao, and iVo = 6000, (b) ajnitial = 0.64ao, 
"collapse ^ — 6.6ao, and TVo = 14500, (c) ajj^j^jg^j = 0.64ao, 
"collapse = — 6.6ao, and A'^o = 5500, and (d) ajnitial = '^"O' 
"collapse ~ — 263ao, and A''o — 6000. Full line: present the- 
ory; dashed line: average over preliminary, unanalyzed data 
using Eq. (0) 



foi" ^evolve > ^collapse- ^^S- 1 we also plot N{t^^^^^^) 
of Eq. ( |3.lD (dashed line) for ^collapse ~ — 263ao, — 30ao 
and — 6.7ao with respective decay rates '''decay ~ ^■'^ 
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2.8 ms and 2.8 ms [g5|. For wide variation of parameters 
"initial' "collapse' ^nd No, r^ecay ^^^^s approximately 
between 1 and 3. The results of the present simulation 
(full line) agree well with the average experimental result 
of Eq. ( |3.lD for three different ^collapse (dashed line) 

Next we repeated our calculation for several other val- 
ues for ajnitiai, '^^collapsc ^^'^ These results are plot- 
ted in Fig. 3 for (a) a-^-^-^^ = 89ao, a^^oUapse = "l^oo, 
and No = 6000, (b) ainitial = 0-64ao, a^oUapse = 
— 6.6ao, and A'o = 14500, (c) ajnitjal = 0.64ao, 
"collapse = -6-6ao, and iVo = 5500, and (d) Oinitial = 
7ao, ^collapse ~ — 263ao, and No = 6000. The agree- 
ment of the result of simulation with unpublished, pre- 
liminary unanalyzed data is good in all four cases re- 
ported in Fig. 3 |Eq]. 




1. 



FIG. 4: The collapse time ^collapse I '^collapse I Z*^" 
'^initial ~ ^' ^'•^ ~ 6000 and ^ = 2. Solid circle with error 
bar: experiment Hh open circle: axially symmetric mean- 
field model of Ref! M ; arrows marked Th and Ex are theo- 
retical (4.49) [L3|| and experimental (3.75) |l(| estimates 
of |icollapse!/^0' respectively, full line: present theory. 



general features in the behavior of remnant number and 
collapse time are discussed in the following. 
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The decay curves in Fig. 3 are different, although they 
have certain general features which determine the de- 
cay constant T^Q^gy, collapse time ^collapse' ^^'^ num- 
ber of atoms in the remnant. Experimentally, the frac- 
tion of atoms that went into the remnant decreased with 
l"collapsel and was ~ 40% for la^oHapsel < 10"o ^^^d 
was ~ 10% for |acollapsel ^ lOOao. Figures 1 and 3 also 
show this behavior. The values of T^QQg^y for plots in 

Figs. 3 (a) — (d) are 1.5 ms, 2.4 ms, 3.3 ms, and 1.9 
ms, respectively, lying in the range ~ 1 — 3 ms p5[ |. The 



FIG. 5: Remnant number vs. initial number for ajj^j^jai = 
7ao and difi'erent Ocollapse ^'^^ (a-) C — 3 and (b) ^ = 2. 
The experimental results with error bars are represented 
by solid triangle, solid circle, solid square, and solid inverted 
triangle for tt^ollapse ~ — 21ao, — 30ao, — lOOao and — 255ao. 
The corresponding theoretical results are represented by open 
triangle, open circle, open square, and open inverted triangle. 

Donley et al. provided a quantitative measurement 
of the variation of collapse time ^collapse ^i^^i the fi- 
nal scattering length Ocollapse ^'-'^ ^ given ctj^j^jj^j^ = 
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and A^o = 6000. We also cal cula ted this variation us- 
ing our model given by Eq. (2.3). In our calculation 
we define ^collapse time at which the spikes (ir- 

regularities), as in Fig. 2, tend to appear in the wave 
function. The results are plotted in Fig. 4 and com- 
pared with experimental data jl^ as well as with an- 
other calculation using the mean-field GP equation in an 
axially symmetric trap The agreement between the 
two theoretical results is very good. There is also qual- 
itative agreement between the experimental data on the 
one hand and the two calculations on the other hand: 
^collapse decreases with |acollapsel/^o starting from an 
infinite value at I a(.Qiig^pgg | — acr for a fixed Nq, that is 
6000 in Fig. 4. For this Nq, ocr is the minimum value of 
I '^collapse I ^^^^ leads to the collapse and explosion. For 
a given Nq, a critical value of n = ncr for collapse can 
be defined via ncr = Noacr/l- As there is discrepancy 
between theoretical and experimental ncr for an axially 
symmetric trap |l^ , the theoretical and experimental 
acr are also supposed to be different. The experimental 
kcT = ncrA^/s = 0.46 and the theoretical kci = 0.55 
[l^ for the axially symmetric trap used in the experi- 
ment at JILA with the asymmetry parameter A = 0.389. 
The theoretical acr should be larger than the experimen- 
tal acr in the same proportion. This might imply that the 
theoretical ^collapse should tend to infinity for a slightly 
larger value of ^collapse ^' 
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FIG. 6: The dimensionless rms sizes a;rms (full line) and yrms 
(dashed line) expressed in units of 1/V2 {I = 26070 A) after 
the jump in the scattering length of a BEG of 16000 ^^Rb 
atoms from oijjjj^jg^j ~ 7ao to ^collapse ~ ~6.7ao as functions 
of time Tg^oive for C = 2. 

Donley et al. measured the number of remnant atoms 
for Ojj^i^jai = 7ao and different initial number Nq and 



'^collapse ^^"^ these results Q are plotted in Figs. 5 (a) 
and (b) and compared with numerical simulation per- 
formed with ^ = 3 and 2, respectively. The agreement is 
good for most cases shown in this figure. For Nq = 6000, 
there is some discrepancy between theoretical and ex- 
perimental remnant numbers. The overall agreement is 
better in the case with ^ = 2 than with f = 3. For 
f = 3 the three-body recombination loss-rate is larger 
and this leads to smaller remnant numbers compared to 
the case with ^ = 2. The theoretical Ncr for a fixed neg- 
ative a^ouapse is given by A^cr = 0.55/A-i/V|acollapsel 
]l2| , 13. For flcoUapse ~ — 255ao, — lOOag, — 30ao, and 
-2Iao, iVcr = 124,317,1057, and 1510, respectively. 
Hence, from Fig. 5 we find that the number in the rem- 
nant could be much larger than A'cr for times on the 
order of tens of milliseconds. However, in our simulation 
such a remnant continues to emit atoms at a much slower 
rate and for very large times on the order of seconds the 
number of atoms eventually tends towards A^cr- 

Donley et al. observed that the remnant condensate 
in all cases oscillated in a highly excited collective state 
with approximate frequencies '^Vq^^IqI and '^i^Yadial being 
predominantly excited. The actual measured frequencies 
are 13.6(6) Hz and 33.4(3) Hz. To find if this behavior 
emerges from the present simulation we plot in Fig. 6 
sizes Xrms and j/rms vs. time for the condensate after 
the jump in the scattering length to — 6.7ao from 7ao for 
Nq = 16000. Excluding the first 20 ms when the remnant 
condensate is being formed, we find a periodic oscillation 
in Xrms and yrms with frequencies 13.5 Hz and 34 Hz, 
respectively, as observed in experiment. 



IV. DISCUSSION 

Though we have explained some aspects of the exper- 
iment at JILA, certain detailed features have not been 
addressed in this study. Donley et al. have classified the 
emitted atoms in three categories: burst, missing (unde- 
tected) and jet atoms |Q. The jet atoms appear with 
much lower energy solely in the radial direction possibly 
from the spikes in the wave function when the collapse is 
suddenly interrupted during the period of atom loss be- 
fore the remnant is formed. Strangely enough the emis- 
sion of jet atoms are found not to possess axial symmetry 
always and hence can not be properly treated in a axi- 
ally symmetric model. Moreover a clear-cut distinction 
between the burst and missing atoms emitted during the 
explosion seems to be difficult in the present model as 
the experiment could not specify the properties (magni- 
tude and direction of velocities) of the missing atoms. 
Also, because of the missing atoms it is difficult to pre- 
dict the energy distribution of the burst atoms during the 
explosion in a mean-field analysis. Without proper iden- 
tification of the missing atoms, any energy distribution 
calculated using the present mean-field analysis will yield 
the total energy of burst plus missing atoms. A careful 
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analysis of the energy of the emitted atoms is required 
for explaining the exclusive features and a detailed study 
of the wave function is needed for this purpose. Such an 
analysis is beyond the scope of the present investigation 
and would be a welcome future theoretical work. 

The success of the Crank-Nicholson algorithm in alter- 
nate directions as used in this study depends on a proper 
discretization of the GP equation in space and time. In 
this study we employed a two-dimensional lattice in space 
of 600 X 150 or 90000 points (a; < 15, -30 < y < 30) 
and a time step of 0.001. In the absence of collapse and 
recombination loss this discretization leads to very pre- 
cise results. The accuracy reduces in the presence of the 
violent collapse and explosion simulated by three-body 
recombination. By varying the space discretization grid 
and time step we found that the estimated error in the 
present calculation is less than ~ 10% for time propaga- 
tion upto few tens of milliseconds. 

There has been another attempt to use the mean-field 
GP equation in an axially symmetric trap |^ to explain 
the experiment of Ref. ||lj|. There are certain differ- 
ences between the analysis of Ref. Q and the present 
investigation. According to the experiment of Ref. , 
the burst atoms and missing atoms are components of 
expelled atoms which lose contact with the central con- 
densate that eventually forms the remnant. Of these the 
burst atoms have energy much less than the magnetic 
trap depth. Hence, though expelled from the central con- 
densate they continue tra ppe d and oscillate with time. 
The wave function of Eq. (2^) only describes the central 
condensate. However, in Ref. [|j the burst atoms are 
considered to be the peripheral part (the spikes) of the 
central condensate and hence taken to be described by 
the mean- field Eq. (2^). The missing atoms are actually 
parts of the expelled atoms that have disappeared from 
the trap Q. In Ref. the missing atoms have been 
taken to be the only component of the emitted atoms. 
These are the main differences between the point of view 
of the present analysis and that of Ref. [|| . 

The three-body loss rates of the two studies are also 
widely different. Here we employ the three-body recom- 
bination loss-rate K3 ~ 9 x 10^^^ cm^/s for a = — 370ao 
whereas in Ref. Q the value K3 ~ 10~^^ cm^/s has been 
considered. The present rate is in rough agreement with 
the experimental rate of Ref. |l|l {K3 - 4.2 x lO^^s 
cm^/s) and with the theoretical rate of Ref. |l^ {K3 ^ 
6.7 X 10~^^ cm^/s) for the same value of scattering length 
whereas that of Ref. is orders of magnitude smaller. 
However, such a small three-body rate in Ref. has 
led to a large residual condensate at large time that they 
have interpreted as the sum of burst plus remnant. The 
use of a large three-body rate in this study has led to a 
much smaller residual central condensate which has been 
identified as the remnant as in the experiment at JILA 
0- 

However, it is assuring to see that the ^collapse 
I ''collapse 1/^0 curve of the two models in Fig. 4 agrees 
with each other. The present calculation in Fig. 4 was 



performed with a nonzero loss rate K^, whereas that in 
Ref. 1^ was performed by setting K3 — 0. We find that 
K3 plays an insignificant role in this calculation at small 
times. Hence, the two computer routines lead to the 
same result in the absence of recombination loss before 
the beginning of the explosion. 



V. CONCLUSION 

In conclusion, we have employed a numerical simula- 
tion based on the accurate solution jl^ of the mean-field 
Gross-Pitaevskii equation with a cylindrical trap to study 
the dynamics of the collapse and explosion as observed in 
the recent experiment at JILA |Q. In the GP equation 
we include a quintic three-body nonlinear recombination 
loss term that accounts for the decay of the strongly at- 
tractive condensate. The results of the present simulation 
accounts for some aspects of the experiment. 

In the experiment a strongly attractive *^Rb conden- 
sate was prepared by ramping the scattering length to 
a large negative value and the subsequent decay of the 
collapsing and exploding condensate was measured. We 
have been able to understand the following features of 
this dynamics from the present numerical simulation: 
(1) The condensate undergoes collapse and explosion 
and finally stabilizes to a remnant condensate contain- 
ing about 10% (for |acollapsel > lOOao ) to 40% (for 
I ''collapse I ^ lO^o) of initial number of atoms iVo at large 
times. This percentage is independent of Nq and the 
ramped scattering length Ocollapse- number in the 
remnant condensate can be much larger than the critical 
number for collapse iVcr for the same atomic interaction 
for experimental times on the order of tens of millisec- 
onds. (2) Both in the experiment and our simulation the 
remnant condensate executes radial and axial oscillations 
in a highly excited collective state for a long time with 
frequencies 2j/j.g^jjg^j and '^Vg^xiaV ('^) ^fter the sudden 
change in the scattering length to a large negative value, 
the condensate needs an interval of time ^collapse before 
it experiences loss via explosion. Consequently, the de- 
cay starts after the interval of time ^collapse' ('^) 
number of atoms in the condensate decays exponentially 
with a decay constant t^^qq^ of few milliseconds (~ 1 — 3 
ms). 

To conclude, a large part of the Bosenova experiment 
on ^^Rb atoms at JILA |l^, specially the detailed be- 
havior of the remnant, can be understood by introducing 
the rather conventional three-body recombination loss in 
the standard mean-field GP equation, with a loss rate 
compatible with other studies ||l^, |l^ . The study of the 
detailed behavior of the burst and missing atoms and the 
formation of the jet in such a mean-field theory seems to 
be more complicated technically, nevertheless viable, and 
would be a subject of future investigation. 
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